#pragma once
#include <cstring>
#include <cmath>
enum unit_type {
  U_LENG, //m
  U_MASS, //kg
  U_TIME, //s
  U_RAD,
  U_NONT
};

struct unit {
  static constexpr double cpow(double value, int n){
    double ret = 1;
    if (n < 0) {
      for (int i = 0; i < -n; i++)
        ret /= value;
    } else {
      for (int i = 0; i < n; i ++)
        ret *= value;
    }
    return ret;
  }
  const double value;
  const int p_leng, p_mass, p_time, p_rad;
  constexpr unit(double value, int p_leng, int p_mass, int p_time, int p_rad) : value(value), p_leng(p_leng), p_mass(p_mass), p_time(p_time), p_rad(p_rad) {
  }
  constexpr unit(double value, int utype) : value(value), p_leng(utype == U_LENG), p_mass(utype == U_MASS), p_time(utype == U_TIME), p_rad(utype == U_RAD) {}
  constexpr unit operator()(int n) const {
    return unit(cpow(value, n), p_leng * n, p_mass * n, p_time * n, p_rad * n);
  }
  constexpr explicit operator double() const {
    return value;
  }
  constexpr bool nontype(){
    return p_leng == 0 && p_mass == 0 && p_time == 0;
  }
};
//   __always_inline constexpr unit operator*(const unit &o) {
//     unit ret;
//     for (int i = 0; i < U_NONT; i ++) {
//       ret.powers[i] = powers[i] + o.powers[i];
//     }
//     ret.value = value * o.value;
//     return ret;
//   }
//   __always_inline constexpr unit operator/(const unit &o) {
//     unit ret;
//     for (int i = 0; i < U_NONT; i ++) {
//       ret.powers[i] = powers[i] - o.powers[i];
//     }
//     ret.value = value / o.value;
//     return ret;
//   }
// };
__always_inline constexpr unit operator*(const unit &a, const unit &o) {
  return unit(a.value * o.value, a.p_leng + o.p_leng, a.p_mass + o.p_mass, a.p_time + o.p_time, a.p_rad + o.p_rad);
}
__always_inline constexpr unit operator/(const unit &a, const unit &o) {
  return unit(a.value / o.value, a.p_leng - o.p_leng, a.p_mass - o.p_mass, a.p_time - o.p_time, a.p_rad - o.p_rad);
}
constexpr __always_inline unit operator*(const unit &u, double val) {
  return u * unit(val, 0, 0, 0, 0);
}
constexpr __always_inline unit operator/(const unit &u, double val) {
  return u / unit(val, 0, 0, 0, 0);
}
constexpr __always_inline unit operator*(double val, const unit &u) {
  return unit(val, 0, 0, 0, 0) * u;
}
constexpr __always_inline unit operator/(double val, const unit &u) {
  return unit(val, 0, 0, 0, 0) / u;
}
// /* Si units */
namespace units{
  constexpr unit one(1, U_NONT);
  constexpr unit m(1.0, U_LENG);
  constexpr unit kg(1.0, U_MASS);
  constexpr unit s(1.0, U_TIME);
  constexpr unit rad(1.0, U_RAD);
  constexpr unit m_s = m * s(-1);
  constexpr unit m_s2 = m * s(-2);
  constexpr unit N = m_s2 * kg;
  constexpr unit J = N * m;
  constexpr unit kJ = 1000 * J;

  constexpr unit Ang = 1e-10 * m;
  constexpr unit nm = 1e-9 * m;
  constexpr unit mole(6.02214076e23, U_NONT);
  constexpr unit amu = kg * 0.001 / mole;
  constexpr unit g_mole = amu;
  constexpr unit kCal = 4184 * J;
  constexpr unit kcal_mole = 4184 * J / mole;
  constexpr unit deg = M_PI / 180 * rad;
  constexpr unit gmx_K(double value, int p_nm, int p_rad, int p_kJ = 1, int p_mole = -1) {
    return value * kJ(p_kJ) * nm(p_nm) * mole(p_mole) * rad(p_rad);
  }
};
constexpr double factor_gmx_K(double value, int p_nm, int p_rad, int p_kJ = 1, int p_mole = -1){
  using namespace units;
    unit gmx_unit = value * kJ(p_kJ) * nm(p_nm) * mole(p_mole) * rad(p_rad);
    unit esmd_unit = kCal(p_kJ) * Ang(p_nm) * mole(p_mole) * rad(p_rad);
    return (esmd_unit / gmx_unit).value;
}
// #include <stdio.h>
// int main(){
//   auto gmx_kb = 0.5*1000*units::J/units::mole/units::nm(2);
//   auto charmm_kb = units::kcal_mole/units::Ang(2);
//   auto gmx_ka = units::gmx_K(0.5, 0, -2);
//   auto charmm_ka = units::kcal_mole/units::rad(2);
//   printf("%g %g\n", gmx_ka.value, charmm_ka.value);
//   printf("%f %d\n", (double)(334.720000 * gmx_ka / charmm_ka), (gmx_ka / charmm_ka).nontype());
//   printf("%d %d %d\n", units::J.p_leng, units::J.p_time, units::J.p_mass);
// }